¶ Statistical
analysis ¶
¶ Agalychnis
lemur ¶
Marcelo Araya-Salas, PhD & Fabiola Chirino ¶
¶ 27-04-2023 ¶
source("~/Dropbox/R_package_testing/brmsish/R/extended_summary.R")
source("~/Dropbox/R_package_testing/brmsish/R/helpers.R")
label_spectro_temp <- function(wave, reference = NULL, detection = NULL, envelope = FALSE,
threshold = NULL, smooth = 5, collevels = seq(-100, 0, 5), palette = viridis::viridis,
template.correlation = NULL, line.x.position = 2, hop.size = NULL, col.line = NULL,
...) {
# adjust wl based on hope.size
if (!is.null(hop.size))
wl <- round(wave@samp.rate * hop.size/1000, 0)
# reset graphic device on exit
oldpar <- par(no.readonly = TRUE)
on.exit(par(oldpar))
if (envelope | !is.null(template.correlation))
par(mfrow = c(2, 1), mar = c(0, 4, 1, 1)) else par(mar = c(4, 4, 1, 1))
# plot spectrogram
seewave::spectro(wave = wave, grid = FALSE, scale = FALSE, palette = palette,
collevels = collevels, axisX = if (envelope | !is.null(template.correlation))
FALSE else TRUE, ...)
# plot detection
if (!is.null(reference))
for (i in seq_len(nrow(reference))) lines(x = (reference[i, c("start", "end")]),
y = rep(line.x.position, 2), col = col.line[1], lwd = 7, lend = 2)
# plot detection
if (!is.null(detection))
for (i in seq_len(nrow(detection))) lines(x = (detection[i, c("start", "end")]),
y = rep(line.x.position - 0.3, 2), col = col.line[2], lwd = 7, lend = 2)
usr <- par("usr")
# add legend
if (!is.null(detection) & !is.null(reference))
legend(x = usr[2] * 0.98, y = usr[4] * 0.98, col = col.line[1:2], legend = c("Referencia",
"Detección"), text.width = 0.3, lwd = 4, bg = "#FFFFFFE6", xjust = 1,
yjust = 1)
if (is.null(detection) & !is.null(reference))
legend(x = usr[2] * 0.98, y = usr[4] * 0.98, col = col.line[1], text.width = 0.3,
legend = c("Referencia"), lwd = 4, bg = "#FFFFFFE6", xjust = 1, yjust = 1)
if (is.null(reference) & !is.null(detection))
legend(x = usr[2] * 0.98, y = usr[4] * 0.98, col = col.line[2], legend = c("Detección"),
lwd = 4, text.width = 0.3, bg = "#FFFFFFE6", xjust = 1, yjust = 1)
if (envelope) {
# set graphic device for envelope
par(mar = c(4, 4, 0.3, 1))
if (!is.null(smooth))
smooth <- round(wave@samp.rate * smooth/1000, 0)
# plot envelope
seewave::env(wave, colwave = "#07889B", ssmooth = smooth)
# add threshold line
if (!is.null(threshold))
abline(h = par("usr")[4] * threshold/100, col = col.line[3], lwd = 3)
} else if (!is.null(template.correlation)) {
# set graphic device for correlations
par(mar = c(4, 4, 0.3, 1))
plot(x = seq(template.correlation$template.duration/2, duration(wave) - template.correlation$template.duration/2,
length.out = length(template.correlation$correlation.scores)), y = template.correlation$correlation.scores,
type = "l", xlab = "Tiempo (s)", ylab = "Correlación", col = "#07889B",
lwd = 1.6, xaxs = "i", xlim = c(0, duration(wave)))
# add threshold line
if (!is.null(threshold))
abline(h = threshold, col = col.line[3], lwd = 3)
}
}clim_dat_2020 <- read_excel("./data/datos_metereologicos/Estacion_met/Dat_met_estacion_2022-01-11.xlsx",
sheet = "2020")
clim_dat_2019 <- read_excel("./data/datos_metereologicos/Estacion_met/Dat_met_estacion_2022-01-11.xlsx",
sheet = "2019")
clim_dat <- rbind(clim_dat_2019, clim_dat_2020)
clim_dat <- clim_dat[, c("filename", "Año", "Mes", "Día", "Hora", "Temp (°C)",
"Humedad Relat.", "Precipitación")]
names(clim_dat) <- c("filename", "year", "month", "day", "hour", "temp", "HR", "rain")
clim_dat <- aggregate(cbind(rain, temp, HR) ~ filename + year + month + day + hour,
clim_dat, mean)
clim_dat$year <- clim_dat$year + 2000
clim_dat$date <- as.Date(paste(clim_dat$year, clim_dat$month, clim_dat$day, sep = "-"))
clim_dat$date_hour <- paste(gsub("-", "", clim_dat$date), clim_dat$hour, sep = "-")
call_dat <- read.csv("./data/processed/call_rate_per_date_time_and_site.csv")
# remove 5 pm data and keep only form Sukia
call_dat <- call_dat[call_dat$hour != 17, ]
call_dat <- call_dat[call_dat$site != "LAGCHIMU", ]
call_dat$date_hour <- paste(sapply(as.character(call_dat$site_date_hour), function(x) strsplit(x,
"_")[[1]][2]), call_dat$hour, sep = "-")
call_dat$temp <- sapply(1:nrow(call_dat), function(x) {
y <- clim_dat$temp[clim_dat$date_hour == call_dat$date_hour[x]]
if (length(y) < 1)
y <- NA
return(y)
})
call_dat$HR <- sapply(1:nrow(call_dat), function(x) {
y <- clim_dat$HR[clim_dat$date_hour == call_dat$date_hour[x]]
if (length(y) < 1)
y <- NA
return(y)
})
call_dat$rain <- sapply(1:nrow(call_dat), function(x) {
y <- clim_dat$rain[clim_dat$date_hour == call_dat$date_hour[x]]
if (length(y) < 1)
y <- NA
return(y)
})
# proportion of acoustic data with climatic data
sum(call_dat$date_hour %in% clim_dat$date_hour)/nrow(call_dat)
sum(!is.na(call_dat$temp))/nrow(call_dat)
call_dat <- call_dat[!is.na(call_dat$temp), ]
call_dat$day <- as.numeric(substr(sapply(as.character(call_dat$site_date_hour), function(x) strsplit(x,
"_")[[1]][2]), 7, 8))
call_dat$date <- as.Date(paste(call_dat$year, call_dat$month, call_dat$day, sep = "-"))
call_dat$moon.date <- ifelse(call_dat$hour < 12, as.Date(call_dat$date - 1), as.Date(call_dat$date))
call_dat$moon.date <- as.Date(call_dat$moon.date, origin = "1970-01-02")
## add moon
call_dat$moonlight <- lunar.illumination(call_dat$moon.date, shift = -6)
call_dat$date_hour_min <- strptime(paste(paste(call_dat$year, call_dat$month, call_dat$day,
sep = "-"), paste(call_dat$hour, "00", sep = ":")), format = "%Y-%m-%d %H:%M")
call_dat$hour_diff <- as.numeric(call_dat$date_hour_min - min(call_dat$date_hour_min))/3600
call_dat$rain_24 <- sapply(1:nrow(call_dat), function(x) sum(clim_dat$rain[strptime(clim_dat$date,
format = "%Y-%m-%d") == (strptime(call_dat$date[x], format = "%Y-%m-%d") - 60 *
60 * 24)]))
call_dat$rain_48 <- sapply(1:nrow(call_dat), function(x) sum(clim_dat$rain[strptime(clim_dat$date,
format = "%Y-%m-%d") == (strptime(call_dat$date[x], format = "%Y-%m-%d") - 60 *
60 * 48)]))
clim_dat$date_hour_min <- strptime(paste(paste(clim_dat$year, clim_dat$month, clim_dat$day,
sep = "-"), paste(clim_dat$hour, "00", sep = ":")), format = "%Y-%m-%d %H:%M")
clim_dat$hour_diff <- as.numeric(clim_dat$date_hour_min - min(call_dat$date_hour_min))/3600
call_dat$prev_temp <- sapply(1:nrow(call_dat), function(x) {
# if(call_dat$hour_diff[x] < 48) pt <- NA else
pt <- mean(clim_dat$temp[clim_dat$hour_diff %in% (call_dat$hour_diff[x] - 48):(call_dat$hour_diff[x] -
24)])
return(pt)
})
write.csv(call_dat, "./data/processed/acoustic_and_climatic_data_by_hour.csv")
Purpose
- Evaluate effect of enviromental factors on vocal activity of A. lemur
Prepare data
Descriptive stats
call_dat_site <- read.csv("./data/processed/call_rate_per_date_time_and_site.csv")- Total number of recordings: 9681
- Total recordings per site:
agg_recs <- aggregate(rec_time ~ site, data = call_dat_site, length)
names(agg_recs)[1:2] <- c("site", "rec_count")
# print table as kable
kb <- kable(agg_recs, row.names = TRUE, digits = 3)
kb <- kable_styling(kb, bootstrap_options = c("striped", "hover", "condensed", "responsive"))
print(kb)| site | rec_count | |
|---|---|---|
| 1 | LAGCHIMU | 5127 |
| 2 | SUKIA | 4554 |
Total recording time: 3224
Total recording time per site:
agg_recs <- aggregate(rec_time ~ site, data = call_dat_site, sum)
names(agg_recs)[1:2] <- c("site", "recording_time")
agg_recs$recording_time <- round((agg_recs$recording_time)/60)
# print table as kable
kb <- kable(agg_recs, row.names = TRUE, digits = 3)
kb <- kable_styling(kb, bootstrap_options = c("striped", "hover", "condensed", "responsive"))
print(kb)| site | recording_time | |
|---|---|---|
| 1 | LAGCHIMU | 1707 |
| 2 | SUKIA | 1517 |
Total detections: 540203
Total detections per site:
agg_recs <- aggregate(n_call ~ site, data = call_dat_site, sum)
names(agg_recs)[1:2] <- c("calls", "count")
# print table as kable
kb <- kable(agg_recs, row.names = TRUE, digits = 3)
kb <- kable_styling(kb, bootstrap_options = c("striped", "hover", "condensed", "responsive"))
print(kb)| calls | count | |
|---|---|---|
| 1 | LAGCHIMU | 169485 |
| 2 | SUKIA | 370718 |
Call rate: 18.598587
Call rate per site:
agg_recs <- aggregate(call_rate ~ site, data = call_dat_site, mean)
agg_recs$sd <- aggregate(call_rate ~ site, data = call_dat_site, sd)[, 2]
names(agg_recs)[1:3] <- c("site", "call_rate", "sd")
# print table as kable
kb <- kable(agg_recs, row.names = TRUE, digits = 3)
kb <- kable_styling(kb, bootstrap_options = c("striped", "hover", "condensed", "responsive"))
print(kb)| site | call_rate | sd | |
|---|---|---|---|
| 1 | LAGCHIMU | 11.018 | 18.039 |
| 2 | SUKIA | 27.133 | 87.778 |
call_rate_hour <- read.csv("./data/processed/acoustic_and_climatic_data_by_hour.csv")
agg <- aggregate(cbind(temp, prev_temp, HR, rain, rain_24, rain_48, moonlight) ~
1, call_rate_hour, function(x) round(c(mean(x), sd(x), min(x), max(x)), 3))
agg <- as.data.frame(matrix(unlist(agg), ncol = 4, byrow = TRUE, dimnames = list(c("Temperature",
"Previous temperature", "Relative humidity", "Night rain", "Rain 24 hours", "Rain 48 hours",
"Moonlight"), c("mean", "sd", "min", "max"))))
# print table as kable
kb <- kable(agg, row.names = TRUE, digits = 3)
kb <- kable_styling(kb, bootstrap_options = c("striped", "hover", "condensed", "responsive"))
print(kb)| mean | sd | min | max | |
|---|---|---|---|---|
| Temperature | 24.229 | 2.126 | 19.472 | 31.763 |
| Previous temperature | 24.097 | 0.988 | 20.005 | 27.111 |
| Relative humidity | 90.056 | 8.464 | 55.158 | 99.940 |
| Night rain | 0.079 | 0.457 | 0.000 | 10.417 |
| Rain 24 hours | 1.441 | 3.161 | 0.000 | 25.620 |
| Rain 48 hours | 1.453 | 3.170 | 0.000 | 25.620 |
| Moonlight | 0.495 | 0.355 | 0.000 | 1.000 |
Mean and sd temperature: 24.229179
Mean previous temperature: 24.097073
Mean temperature: 24.229179
Mean cumulative rain per hour: 0.079232
Mean cumulative rain per hour previous 24 hours: 1.44064
Mean daily cumulative rain per hour previous 48 hours: 1.45317
Bayesian regression models
Scale variables and set model parameters
call_rate_hour <- read.csv("./data/processed/acoustic_and_climatic_data_by_hour.csv")
# make hour a factor
call_rate_hour$hour <- factor(call_rate_hour$hour)
# scale and mean-center
call_rate_hour$sc_temp <- scale(call_rate_hour$temp)
call_rate_hour$sc_HR <- scale(call_rate_hour$HR)
call_rate_hour$sc_rain <- scale(call_rate_hour$rain)
call_rate_hour$sc_rain_24 <- scale(call_rate_hour$rain_24)
call_rate_hour$sc_rain_48 <- scale(call_rate_hour$rain_48)
call_rate_hour$sc_moonlight <- scale(call_rate_hour$moonlight)
call_rate_hour$sc_prev_temp <- scale(call_rate_hour$prev_temp)
priors <- c(prior(normal(0, 4), class = "b"))
chains <- 4
iter <- 10000Global models
The increase in relative humidity, decrease in temperature, increase in the previous accumulated rain, decrease in the night rain, decrease in the percentage of the moon illuminated cause the activity of A. lemur to increase?
fit_glob1 <- brm(n_call | resp_rate(rec_time) ~ sc_temp + sc_HR + sc_moonlight +
sc_rain + sc_rain_24 + sc_prev_temp + ar(p = 2, time = hour_diff, gr = hour),
data = call_rate_hour, iter = iter, chains = chains, cores = chains, family = negbinomial(),
prior = priors, file = "./data/processed/regression_models/global_rain_24", file_refit = "always")
fit_glob2 <- brm(n_call | resp_rate(rec_time) ~ sc_temp + sc_HR + sc_moonlight +
sc_rain + sc_rain_48 + sc_prev_temp + ar(p = 2, time = hour_diff, gr = hour),
data = call_rate_hour, iter = iter, chains = chains, cores = chains, family = negbinomial(),
prior = priors, file = "./data/processed/regression_models/global_rain_48", file_refit = "always")extended_summary(read.file = "./data/processed/regression_models/global_rain_24.rds",
n.posterior = 1000)global_rain_24
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_temp + sc_HR + sc_moonlight + sc_rain + sc_rain_24 + sc_prev_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 | 0 | 9377.5 | 12503.8 | 536473866 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_Intercept | 0.698 | 0.636 | 0.760 | 1 | 9377.5 | 12503.8 |
| b_sc_temp | 0.573 | 0.493 | 0.654 | 1 | 10930.3 | 12951.7 |
| b_sc_HR | 0.307 | 0.233 | 0.381 | 1 | 14333.7 | 14034.8 |
| b_sc_moonlight | -0.090 | -0.149 | -0.031 | 1 | 12403.9 | 13102.6 |
| b_sc_rain | 0.038 | 0.002 | 0.075 | 1 | 18321.1 | 15080.0 |
| b_sc_rain_24 | 0.091 | 0.053 | 0.130 | 1 | 18798.2 | 14951.0 |
| b_sc_prev_temp | -0.024 | -0.069 | 0.022 | 1 | 14546.5 | 13679.8 |
extended_summary(read.file = "./data/processed/regression_models/global_rain_48.rds",
n.posterior = 1000)global_rain_48
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_temp + sc_HR + sc_moonlight + sc_rain + sc_rain_48 + sc_prev_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 | 0 | 12375.1 | 13473.4 | 58053599 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_Intercept | 0.703 | 0.640 | 0.766 | 1 | 12375.1 | 13699.9 |
| b_sc_temp | 0.566 | 0.485 | 0.649 | 1 | 13593.6 | 13473.4 |
| b_sc_HR | 0.303 | 0.228 | 0.379 | 1 | 17757.9 | 15101.0 |
| b_sc_moonlight | -0.098 | -0.158 | -0.039 | 1 | 18857.2 | 14680.5 |
| b_sc_rain | 0.031 | -0.006 | 0.068 | 1 | 22865.9 | 15592.2 |
| b_sc_rain_48 | -0.006 | -0.044 | 0.033 | 1 | 24427.4 | 16403.0 |
| b_sc_prev_temp | -0.037 | -0.084 | 0.010 | 1 | 20592.6 | 15806.9 |
Temperature
The increase in temperature at night causes that activity of A. lemur to decrease?
fit2.1 <- brm(n_call | resp_rate(rec_time) ~ sc_temp + sc_rain + sc_rain_24 + ar(p = 2,
time = hour_diff, gr = hour), data = call_rate_hour, iter = iter, chains = chains,
cores = chains, family = negbinomial(), prior = priors, file = "./data/processed/regression_models/temp2.1",
file_refit = "always")
fit2.2 <- brm(n_call | resp_rate(rec_time) ~ sc_temp + sc_rain + sc_rain_48 + ar(p = 2,
time = hour_diff, gr = hour), data = call_rate_hour, iter = iter, chains = chains,
cores = chains, family = negbinomial(), prior = priors, file = "./data/processed/regression_models/temp2.2",
file_refit = "always")
fit2.3 <- brm(n_call | resp_rate(rec_time) ~ sc_prev_temp + ar(p = 2, time = hour_diff,
gr = hour), data = call_rate_hour, iter = iter, chains = chains, cores = chains,
family = negbinomial(), prior = priors, file = "./data/processed/regression_models/temp2.3",
file_refit = "always")extended_summary(read.file = "./data/processed/regression_models/temp2.1.rds", n.posterior = 1000)temp2.1
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_temp + sc_rain + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 1 | 0 | 10328.4 | 12111.8 | 1846130948 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_Intercept | 0.698 | 0.634 | 0.762 | 1 | 10328.4 | 12111.8 |
| b_sc_temp | 0.312 | 0.261 | 0.361 | 1 | 12331.5 | 14350.8 |
| b_sc_rain | 0.044 | 0.008 | 0.083 | 1 | 20658.2 | 15981.9 |
| b_sc_rain_24 | 0.090 | 0.052 | 0.128 | 1 | 20542.9 | 16470.9 |
extended_summary(read.file = "./data/processed/regression_models/temp2.2.rds", n.posterior = 1000)temp2.2
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_temp + sc_rain + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 | 0 | 8873.8 | 11982.2 | 2046380001 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_Intercept | 0.701 | 0.635 | 0.766 | 1 | 8873.8 | 11982.2 |
| b_sc_temp | 0.305 | 0.255 | 0.355 | 1 | 10391.0 | 13564.9 |
| b_sc_rain | 0.036 | 0.000 | 0.073 | 1 | 17117.4 | 15242.3 |
| b_sc_rain_48 | 0.006 | -0.031 | 0.044 | 1 | 17732.8 | 15884.1 |
extended_summary(read.file = "./data/processed/regression_models/temp2.3.rds", n.posterior = 1000)temp2.3
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_prev_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 | 0 | 20200 | 15620.2 | 1962664000 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_Intercept | 0.675 | 0.605 | 0.745 | 1 | 20200.0 | 15620.2 |
| b_sc_prev_temp | -0.002 | -0.048 | 0.044 | 1 | 36045.7 | 16766.2 |
Relative humidity
Does an increase in relative humidity cause the activity of A. lemur to increase?
fit.3.1 <- brm(n_call | resp_rate(rec_time) ~ sc_temp + sc_HR + sc_rain + sc_rain_24 +
ar(p = 2, time = hour_diff, gr = hour), data = call_rate_hour, iter = iter, chains = chains,
cores = chains, family = negbinomial(), prior = priors, file = "./data/processed/regression_models/RH3.1",
file_refit = "always")
fit.3.2 <- brm(n_call | resp_rate(rec_time) ~ sc_temp + sc_HR + sc_rain + sc_rain_48 +
ar(p = 2, time = hour_diff, gr = hour), data = call_rate_hour, iter = iter, chains = chains,
cores = chains, family = negbinomial(), prior = priors, file = "./data/processed/regression_models/RH3.2",
file_refit = "always") # Este tengo que doblemente confirmarloextended_summary(read.file = "./data/processed/regression_models/RH3.1.rds", n.posterior = 1000)RH3.1
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_temp + sc_HR + sc_rain + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 | 0 | 10000 | 12606.5 | 1854758574 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_Intercept | 0.699 | 0.637 | 0.760 | 1 | 10000.0 | 12606.5 |
| b_sc_temp | 0.578 | 0.499 | 0.656 | 1 | 13756.3 | 13884.5 |
| b_sc_HR | 0.314 | 0.242 | 0.385 | 1 | 17498.1 | 14651.1 |
| b_sc_rain | 0.037 | 0.002 | 0.075 | 1 | 24287.1 | 16414.2 |
| b_sc_rain_24 | 0.096 | 0.058 | 0.135 | 1 | 22809.6 | 15579.0 |
extended_summary(read.file = "./data/processed/regression_models/RH3.2.rds", n.posterior = 1000)RH3.2
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_temp + sc_HR + sc_rain + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 | 0 | 8492.95 | 12550.7 | 1444793977 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_Intercept | 0.703 | 0.641 | 0.765 | 1 | 8492.95 | 12697.7 |
| b_sc_temp | 0.567 | 0.488 | 0.648 | 1 | 9626.63 | 12550.7 |
| b_sc_HR | 0.308 | 0.235 | 0.382 | 1 | 12079.96 | 13794.6 |
| b_sc_rain | 0.029 | -0.007 | 0.066 | 1 | 18138.82 | 15107.9 |
| b_sc_rain_48 | 0.005 | -0.032 | 0.043 | 1 | 16926.99 | 14915.2 |
Moon
Decreasing the percentage of the moon illuminated causes an increase in A. lemur activity?
fit.4 <- brm(n_call | resp_rate(rec_time) ~ sc_moonlight + ar(p = 2, time = hour_diff,
gr = hour), data = call_rate_hour, iter = iter, chains = chains, cores = chains,
family = negbinomial(), prior = priors, file = "./data/processed/regression_models/moon4",
file_refit = "always")extended_summary(read.file = "./data/processed/regression_models/moon4.rds", n.posterior = 1000)moon4
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 | 0 | 19357 | 14632 | 1734290538 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_Intercept | 0.674 | 0.606 | 0.743 | 1 | 19357.0 | 14632.0 |
| b_sc_moonlight | -0.128 | -0.194 | -0.063 | 1 | 31443.5 | 15209.9 |
Night rain
If the night rain decreases can affect the A. lemur activity to increase?
fit.5.1 <- brm(n_call | resp_rate(rec_time) ~ sc_prev_temp + sc_rain + sc_rain_48 +
ar(p = 2, time = hour_diff, gr = hour), data = call_rate_hour, iter = iter, chains = chains,
cores = chains, family = negbinomial(), prior = priors, file = "./data/processed/regression_models/night_rain5.1",
file_refit = "always")
fit.5.2 <- brm(n_call | resp_rate(rec_time) ~ sc_prev_temp + sc_rain + sc_rain_24 +
ar(p = 2, time = hour_diff, gr = hour), data = call_rate_hour, iter = iter, chains = chains,
cores = chains, family = negbinomial(), prior = priors, file = "./data/processed/regression_models/night_rain5.2",
file_refit = "always")extended_summary(read.file = "./data/processed/regression_models/night_rain5.1.rds",
n.posterior = 1000)night_rain5.1
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_prev_temp + sc_rain + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 | 0 | 15449.5 | 13446.5 | 1903747053 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_Intercept | 0.676 | 0.607 | 0.745 | 1.001 | 15449.5 | 13446.5 |
| b_sc_prev_temp | 0.003 | -0.046 | 0.052 | 1 | 27583.8 | 14753.3 |
| b_sc_rain | -0.006 | -0.042 | 0.030 | 1 | 29045.4 | 16584.0 |
| b_sc_rain_48 | 0.014 | -0.027 | 0.054 | 1 | 27750.2 | 15632.5 |
extended_summary(read.file = "./data/processed/regression_models/night_rain5.2.rds",
n.posterior = 1000)night_rain5.2
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_prev_temp + sc_rain + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 | 0 | 21570 | 14952.6 | 1800703315 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_Intercept | 0.673 | 0.604 | 0.742 | 1 | 21570.0 | 14952.6 |
| b_sc_prev_temp | 0.007 | -0.039 | 0.054 | 1 | 39064.3 | 16543.7 |
| b_sc_rain | 0.001 | -0.035 | 0.038 | 1 | 36294.6 | 16198.6 |
| b_sc_rain_24 | 0.077 | 0.038 | 0.117 | 1 | 34404.5 | 17412.9 |
Previous Rain
Does an increase in the accumulated previous rain causes that A. lemur activity to increase?
fit.6.1 <- brm(n_call | resp_rate(rec_time) ~ sc_temp + sc_rain + sc_rain_48 + ar(p = 2,
time = hour_diff, gr = hour), data = call_rate_hour, iter = iter, chains = chains,
cores = chains, family = negbinomial(), prior = priors, file = "./data/processed/regression_models/previous_rain6.1",
file_refit = "always")
fit.6.2 <- brm(n_call | resp_rate(rec_time) ~ sc_temp + sc_rain + sc_rain_24 + ar(p = 2,
time = hour_diff, gr = hour), data = call_rate_hour, iter = iter, chains = chains,
cores = chains, family = negbinomial(), prior = priors, file = "./data/processed/regression_models/previous_rain6.2",
file_refit = "always")
fit.6.3 <- brm(n_call | resp_rate(rec_time) ~ sc_temp + sc_HR + sc_rain_24 + ar(p = 2,
time = hour_diff, gr = hour), data = call_rate_hour, iter = iter, chains = chains,
cores = chains, family = negbinomial(), prior = priors, file = "./data/processed/regression_models/night_rain6.3",
file_refit = "always")extended_summary(read.file = "./data/processed/regression_models/previous_rain6.1.rds",
n.posterior = 1000)previous_rain6.1
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_temp + sc_rain + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 | 0 | 8446.05 | 11423.8 | 20512253 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_Intercept | 0.701 | 0.636 | 0.764 | 1 | 8446.05 | 11423.8 |
| b_sc_temp | 0.304 | 0.254 | 0.355 | 1 | 10470.72 | 13123.6 |
| b_sc_rain | 0.036 | 0.000 | 0.074 | 1 | 16919.56 | 15417.1 |
| b_sc_rain_48 | 0.007 | -0.031 | 0.045 | 1 | 17069.15 | 15935.4 |
extended_summary(read.file = "./data/processed/regression_models/previous_rain6.2.rds",
n.posterior = 1000)previous_rain6.2
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_temp + sc_rain + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 | 0 | 11468.2 | 13337.9 | 1714781937 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_Intercept | 0.698 | 0.633 | 0.762 | 1 | 11468.2 | 13337.9 |
| b_sc_temp | 0.312 | 0.262 | 0.362 | 1 | 13834.2 | 13973.2 |
| b_sc_rain | 0.045 | 0.008 | 0.083 | 1 | 23413.4 | 15796.1 |
| b_sc_rain_24 | 0.090 | 0.051 | 0.129 | 1 | 22867.5 | 15487.7 |
extended_summary(read.file = "./data/processed/regression_models/night_rain6.3.rds",
n.posterior = 1000)night_rain6.3
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_temp + sc_HR + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 2 | 0 | 7697.31 | 11401.3 | 2086406674 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_Intercept | 0.699 | 0.636 | 0.760 | 1 | 7697.31 | 11401.3 |
| b_sc_temp | 0.572 | 0.493 | 0.650 | 1 | 8657.91 | 12183.9 |
| b_sc_HR | 0.318 | 0.244 | 0.391 | 1 | 10071.18 | 12501.7 |
| b_sc_rain_24 | 0.093 | 0.055 | 0.131 | 1 | 14992.10 | 14450.6 |
Causal model combined effect size plot
Takes the effect sizes (and posteriors) from the right causal models
coef_table <- read.csv("./data/processed/summary_causal_model_table_13-01-23.csv",
sep = "\t")
coef_table$variable <- coef_table$Label
coef_table$value <- coef_table$Estimate
coef_table$significance <- ifelse(coef_table$CI_low * coef_table$CI_high > 0, "sig",
"no.sig")
posteriors_l <- lapply(1:nrow(coef_table), function(x) {
# print(x)
X <- readRDS(coef_table$Model[x])
xdrws <- brms::as_draws(X)
post <- xdrws$`1`[[paste0("b_", coef_table$Variable[x])]]
out <- data.frame(variable = coef_table$Label[x], value = post, significance = coef_table$significance[x])
return(out)
})
posteriors <- do.call(rbind, posteriors_l)
coef_table$variable <- factor(coef_table$Label, levels = sort(unique(coef_table$Label),
FALSE))
posteriors$variable <- factor(posteriors$variable, levels = sort(unique(posteriors$variable),
TRUE))
fill_values <- c("#FFB9DF", "#FF7598")
fill_values <- adjustcolor(fill_values, alpha.f = 0.5)
# creat plots gg_dists <-
ggplot2::ggplot(data = posteriors, ggplot2::aes(y = variable, x = value, fill = significance)) +
ggplot2::geom_vline(xintercept = 0, col = "black", lty = 2) + ggdist::stat_halfeye(ggplot2::aes(x = value),
.width = c(0.95), normalize = "panels", color = "transparent") + ggplot2::scale_fill_manual(values = fill_values,
guide = "none") + ggplot2::geom_point(data = coef_table) + ggplot2::geom_errorbar(data = coef_table,
ggplot2::aes(xmin = CI_low, xmax = CI_high), width = 0) + ggplot2::scale_color_manual(values = coef) +
ggplot2::facet_wrap(~variable, scales = "free_y", ncol = 1) + ggplot2::theme_classic() +
ggplot2::theme(axis.ticks.length = ggplot2::unit(0, "pt"), plot.margin = ggplot2::margin(0,
0, 0, 0, "pt"), legend.position = "none", strip.background = ggplot2::element_blank(),
strip.text = ggplot2::element_blank()) + ggplot2::labs(x = "Effect size",
y = "Parameter") #+
# ggplot2::xlim(range(c(posteriors_by_chain$value, 0)) * plot.area.prop)
ggsave(filename = "./figures/summary_effect_sizes_pooled_from_multiple_causal_models.jpeg")Conditional plots
Measured based on the global model
Rain and temperature
glob_mod_24 <- readRDS("./data/processed/regression_models/global_rain_24.rds")
conditions <- data.frame(sc_temp = c(`Low temperature` = -1, `Mean temperature` = 0,
`High temperature` = 1))
rain_gg <- plot(conditional_effects(glob_mod_24, effects = "sc_rain", conditions = conditions),
plot = FALSE)[[1]] + scale_color_viridis_d() + theme_classic() + ylim(c(0, 520)) +
ggtitle("Current rain") + labs(x = "Rain", y = "Call activity (calls/hour)")
rain24_gg <- plot(conditional_effects(glob_mod_24, effects = "sc_rain_24", conditions = conditions),
plot = FALSE)[[1]] + scale_color_viridis_d() + theme_classic() + ylim(c(0, 520)) +
ggtitle("Previous 24h rain") + labs(x = "Rain", y = "Call activity (calls/hour)")
glob_mod_48 <- readRDS("./data/processed/regression_models/global_rain_48.rds")
rain48_gg <- plot(conditional_effects(glob_mod_48, effects = "sc_rain_48", conditions = conditions),
plot = FALSE)[[1]] + scale_color_viridis_d() + theme_classic() + ylim(c(0, 520)) +
ggtitle("Previous 48h rain") + labs(x = "Rain", y = "Call activity (calls/hour)")
cowplot::plot_grid(rain_gg, rain24_gg, rain48_gg, nrow = 3)
Relative humidity and temperature
conditions <- data.frame(sc_temp = c(`Low temperature` = -1, `Mean temperature` = 0,
`High temperature` = 1))
hr_temp_gg <- plot(conditional_effects(glob_mod_24, effects = "sc_HR", conditions = conditions),
plot = FALSE)[[1]] + scale_color_viridis_d() + theme_classic() + labs(x = "Relative humidity",
y = "Call activity (calls/hour)")
hr_temp_gg
conditions <- data.frame(sc_HR = c(`Low humidity` = -1, `Mean humidity` = 0, `High humidity` = 1))
temp_hr_gg <- plot(conditional_effects(glob_mod_24, effects = "sc_temp", conditions = conditions),
plot = FALSE)[[1]] + scale_color_viridis_d() + theme_classic() + labs(x = "Temperature",
y = "Call activity (calls/hour)")
temp_hr_gg
Relative humidity and rain
conditions <- data.frame(sc_rain = c(`Low rain` = -1, `Mean rain` = 0, `High rain` = 1))
rain_gg <- plot(conditional_effects(glob_mod_24, effects = "sc_HR", conditions = conditions),
plot = FALSE)[[1]] + scale_color_viridis_d() + theme_classic() + ylim(c(0, 75)) +
ggtitle("Current rain")
rain24_gg <- plot(conditional_effects(glob_mod_24, effects = "sc_HR", conditions = conditions),
plot = FALSE)[[1]] + scale_color_viridis_d() + theme_classic() + ylim(c(0, 75)) +
ggtitle("Previous 24h rain")
rain48_gg <- plot(conditional_effects(glob_mod_48, effects = "sc_HR", conditions = conditions),
plot = FALSE)[[1]] + scale_color_viridis_d() + theme_classic() + ylim(c(0, 75)) +
ggtitle("Previous 48h rain")
cowplot::plot_grid(rain_gg, rain24_gg, rain48_gg, nrow = 3)
conditions <- data.frame(sc_HR = c(`Low humidity` = -1, `Mean humidity` = 0, `High humidity` = 1))
rain_hr_gg <- plot(conditional_effects(glob_mod_24, effects = "sc_rain", conditions = conditions),
plot = FALSE)[[1]] + scale_color_viridis_d() + theme_classic() + ggtitle("Current rain") +
labs(x = "Rain", y = "Call activity (calls/hour)")
rain24_hr_gg <- plot(conditional_effects(glob_mod_24, effects = "sc_rain_24", conditions = conditions),
plot = FALSE)[[1]] + scale_color_viridis_d() + theme_classic() + ggtitle("Previous 24h rain") +
labs(x = "Rain previous", y = "Call activity (calls/hour)")
rain48_hr_gg <- plot(conditional_effects(glob_mod_48, effects = "sc_rain_48", conditions = conditions),
plot = FALSE)[[1]] + scale_color_viridis_d() + theme_classic() + ggtitle("Previous 48h rain") +
labs(x = "Rain previous", y = "Call activity (calls/hour)")
cowplot::plot_grid(rain_hr_gg, rain24_hr_gg, rain48_hr_gg, nrow = 3)
Moonlight and temperature
conditions <- data.frame(sc_temp = c(`Low temperature` = -1, `Mean temperature` = 0,
`High temperature` = 1))
moon_temp_gg <- plot(conditional_effects(glob_mod_24, effects = "sc_moonlight", conditions = conditions),
plot = FALSE)[[1]] + scale_color_viridis_d() + theme_classic() + labs(x = "Moonlight",
y = "Call activity (calls/hour)")
moon_temp_gg
conditions <- data.frame(sc_moonlight = c(`Low moonlight` = -1, `Mean moonlight` = 0,
`High moonlight` = 1))
temp_moon_gg <- plot(conditional_effects(glob_mod_24, effects = "sc_temp", conditions = conditions),
plot = FALSE)[[1]] + scale_color_viridis_d() + theme_classic() + labs(x = "Temperature",
y = "Call activity (calls/hour)")
temp_moon_gg
# st <- data.frame(sound.files = 'FIGURA_LEMUR.wav', selec = 1, start = 0.1,
# end = 0.4)
# tweak_spectro(st, length.out = 20, ovlp = 99, wl = c(100, 1000), pal =
# c('reverse.gray.colors.2', 'viridis', 'reverse.terrain.colors'), path =
# './figures/Figura_canto_lemur', flim = c(1, 4), ncol = 10, nrow = 6, width =
# 15, collev.min = c(-110))
wav <- readWave("./figures/Figura_canto_lemur/FIGURA_LEMUR.wav")
graphics.off()
par(mfrow = c(2, 1), mar = c(0, 4, 4, 1))
warbleR:::spectro_wrblr_int2(wav, flim = c(1, 4), wl = 700, grid = F, collevels = seq(-125,
0, 1), ovlp = 0, palette = viridis, axisX = FALSE)
peaks <- c(0.320165, 0.9, 1.7, 2.6, 3.5, 4.2, 4.5) # + seq(-0.1, 0.2, length.out = 7)
valleys <- seq(0, duration(wav) + 0.3, length.out = 100)
cor_dat <- data.frame(time = c(peaks, valleys), cor = c(rep(0.7, length(peaks)),
rep(0.1, length(valleys))))
cor_dat$cor <- cor_dat$cor + rnorm(nrow(cor_dat), sd = 0.03)
cor_dat$cor <- smoothw(cor_dat$cor, wl = 2, f = 10)
cor_dat <- cor_dat[order(cor_dat$time), ]
par(mar = c(2, 4, 0, 1))
plot(cor_dat$time, cor_dat$cor, type = "l", xaxs = "i")
spectro(wav, flim = c(1, 4), wl = 700, grid = F, collevels = seq(-125, 0, 1), ovlp = 99,
palette = viridis, axisX = FALSE, tlim = c(1.53, 1.63), scale = FALSE, flab = "Frecuencia (kHz)",
tlab = "Tiempo (s)")template <- data.frame(sound.files = "FIGURA_LEMUR.wav", selec = 1, start = 1.55,
end = 1.61, bottom.freq = 2, top.freq = 3)
# get correlations
correlations <- template_correlator(templates = template, files = "FIGURA_LEMUR.wav",
path = "./figures/Figura_canto_lemur/")
thresh <- 0.7
# run detection
detection <- template_detector(template.correlations = correlations, threshold = thresh)
reference <- template_detector(template.correlations = correlations, threshold = 0.55)
detection
# plot spectrogram
label_spectro_temp(wave = wav, reference = reference, detection = detection, template.correlation = correlations[[1]],
flim = c(1, 4), threshold = thresh, hop.size = 10, ovlp = 50, collevels = seq(-125,
0, 1), col.line = c("#AFBF35", "#F20505", "#F2622E"), flab = "Frecuencia (kHz)",
tlab = "Tiempo (s)")
# 012623 034941 AFBF35 F2622E F20505
Takeaways
Sum up results
Next steps
Session information
## R version 4.1.0 (2021-05-18)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
## LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3
##
## locale:
## [1] LC_CTYPE=pt_BR.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=es_CR.UTF-8 LC_COLLATE=pt_BR.UTF-8
## [5] LC_MONETARY=es_CR.UTF-8 LC_MESSAGES=pt_BR.UTF-8
## [7] LC_PAPER=es_CR.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ohun_0.1.0 warbleR_1.1.28 NatureSounds_1.0.4 seewave_2.2.0
## [5] tuneR_1.4.1 brmsish_1.0.0 lunar_0.2-1 ggplot2_3.4.0
## [9] knitr_1.42 kableExtra_1.3.4 HDInterval_0.2.2 readxl_1.3.1
## [13] posterior_1.3.1 cowplot_1.1.1 ggdist_3.2.0 brms_2.18.0
## [17] Rcpp_1.0.9 viridis_0.6.2 viridisLite_0.4.1 remotes_2.4.2
##
## loaded via a namespace (and not attached):
## [1] backports_1.4.1 systemfonts_1.0.4 plyr_1.8.7
## [4] igraph_1.3.5 splines_4.1.0 crosstalk_1.2.0
## [7] TH.data_1.1-0 rstantools_2.2.0 inline_0.3.19
## [10] digest_0.6.31 htmltools_0.5.4 fansi_1.0.3
## [13] magrittr_2.0.3 checkmate_2.1.0 RcppParallel_5.1.5
## [16] matrixStats_0.62.0 xts_0.12.2 sandwich_3.0-1
## [19] svglite_2.1.0 prettyunits_1.1.1 colorspace_2.0-3
## [22] signal_0.7-7 rvest_1.0.3 textshaping_0.3.5
## [25] xfun_0.36 dplyr_1.0.10 RCurl_1.98-1.9
## [28] callr_3.7.3 crayon_1.5.2 jsonlite_1.8.4
## [31] lme4_1.1-27.1 survival_3.2-11 zoo_1.8-11
## [34] ape_5.6-2 glue_1.6.2 gtable_0.3.1
## [37] emmeans_1.8.1-1 webshot_0.5.4 distributional_0.3.1
## [40] pkgbuild_1.4.0 rstan_2.21.7 abind_1.4-5
## [43] scales_1.2.1 mvtnorm_1.1-3 DBI_1.1.1
## [46] miniUI_0.1.1.1 dtw_1.23-1 xtable_1.8-4
## [49] diffobj_0.3.4 proxy_0.4-27 stats4_4.1.0
## [52] StanHeaders_2.21.0-7 DT_0.26 htmlwidgets_1.5.4
## [55] httr_1.4.4 threejs_0.3.3 ellipsis_0.3.2
## [58] pkgconfig_2.0.3 loo_2.4.1.9000 farver_2.1.1
## [61] sass_0.4.5 utf8_1.2.2 labeling_0.4.2
## [64] tidyselect_1.2.0 rlang_1.0.6 reshape2_1.4.4
## [67] later_1.3.0 munsell_0.5.0 cellranger_1.1.0
## [70] tools_4.1.0 cachem_1.0.6 cli_3.6.0
## [73] generics_0.1.3 ggridges_0.5.4 evaluate_0.20
## [76] stringr_1.5.0 fastmap_1.1.0 ragg_1.1.3
## [79] oce_1.7-8 yaml_2.3.7 processx_3.8.0
## [82] pbapply_1.6-0 nlme_3.1-152 mime_0.12
## [85] projpred_2.0.2 formatR_1.11 rstanarm_2.21.3
## [88] xml2_1.3.3 compiler_4.1.0 bayesplot_1.9.0
## [91] shinythemes_1.2.0 rstudioapi_0.14 gamm4_0.2-6
## [94] tibble_3.1.8 bslib_0.4.2 stringi_1.7.12
## [97] highr_0.10 ps_1.7.2 Brobdingnag_1.2-9
## [100] lattice_0.20-44 Matrix_1.5-1 nloptr_1.2.2.2
## [103] markdown_1.3 fftw_1.0-7 shinyjs_2.1.0
## [106] tensorA_0.36.2 vctrs_0.5.2 pillar_1.8.1
## [109] lifecycle_1.0.3 gsw_1.0-6 jquerylib_0.1.4
## [112] bridgesampling_1.1-2 estimability_1.4.1 bitops_1.0-7
## [115] Sim.DiffProc_4.8 httpuv_1.6.6 R6_2.5.1
## [118] promises_1.2.0.1 gridExtra_2.3 codetools_0.2-18
## [121] boot_1.3-28 colourpicker_1.2.0 MASS_7.3-54
## [124] gtools_3.9.3 assertthat_0.2.1 rjson_0.2.21
## [127] withr_2.5.0 Deriv_4.1.3 shinystan_2.6.0
## [130] multcomp_1.4-17 mgcv_1.8-36 parallel_4.1.0
## [133] grid_4.1.0 coda_0.19-4 minqa_1.2.4
## [136] rmarkdown_2.20 shiny_1.7.3 base64enc_0.1-3
## [139] dygraphs_1.1.1.6